# name : BE (HVD) analysis using proc glm with 3-way design 3x3
# key : be.hvd.glm.3way.
# contributor: Shuguang Sun
# --
* Dataset containing TEST observations;
data test;
set pk;
if trt = 'T';
latt = lauct;
run;

* Dataset containing REFERENCE 1 observations;
data ref1;
set pk;
if trt = 'R';
if (seq = 1 and per = 2) or (seq = 2 and per = 1) or (seq = 3 and per = 1);
lat1r = lauct;
run;

* Dataset containing REFERENCE 2 observations;
data ref2;
set pk;
if trt = 'R';
if (seq = 1 and per = 3) or (seq = 2 and per = 3) or (seq = 3 and per = 2);
lat2r = lauct;
run;

* Determine Iij and Dij;
data scavbe;
merge test ref1 ref2;
by seq subj;
ilat = latt - (0.5 * (lat1r + lat2r));
dlat = lat1r - lat2r;
run;

* Intermediate analysis -- ilat;
proc glm data=scavbe;
class seq;
model ilat = seq / clparm alpha=0.1;
estimate 'average' intercept 1 seq 0.3333333333 0.3333333333 0.3333333333;
ods output overallanova=iglm1;
ods output Estimates=iglm2;
ods output NObs=iglm3;
title1 'scaled average BE';
run;

* From the dataset IGLM2, calculate the following IGLM2;
data iglm2;
set iglm2;
pointest = exp(estimate);
x = estimate**2 - stderr**2;
boundx = (max((abs(LowerCL)), (abs(UpperCL))))**2;
run;

* Intermediate analysis -- dlat;
proc glm data=scavbe;
class seq;
model dlat = seq;
ods output overallanova=dglm1;
ods output NObs=dglm3;
title1 'scaled average BE';
run;

* From the dataset DGLM1, calculate the following DGLM1;
data dglm1;
set dglm1;
dfd = df;
s2wr = ms / 2;
run;

* From the above parameters, calculate the final 95 % upper confidence bound;
data result;
merge iglm2 dglm1;
theta = ((log(1.25)) / 0.25)**2;
y = - theta * s2wr;
boundy = y * dfd / cinv(0.95, dfd);
sWR = sqrt(s2wr);
critbound = (x + y) + sqrt(((boundx - x)**2) + ((boundy - y)**2));
run;